//-------------------------------------------------------------------------//
//                                                                         //
//  This benchmark is an OpenMP C version of the NPB SP code. This OpenMP  //
//  C version is developed by the Center for Manycore Programming at Seoul //
//  National University and derived from the OpenMP Fortran versions in    //
//  "NPB3.3-OMP" developed by NAS.                                         //
//                                                                         //
//  Permission to use, copy, distribute and modify this software for any   //
//  purpose with or without fee is hereby granted. This software is        //
//  provided "as is" without express or implied warranty.                  //
//                                                                         //
//  Information on NPB 3.3, including the technical report, the original   //
//  specifications, source code, results and information on how to submit  //
//  new results, is available at:                                          //
//                                                                         //
//           http://www.nas.nasa.gov/Software/NPB/                         //
//                                                                         //
//  Send comments or suggestions for this OpenMP C version to              //
//  cmp@aces.snu.ac.kr                                                     //
//                                                                         //
//          Center for Manycore Programming                                //
//          School of Computer Science and Engineering                     //
//          Seoul National University                                      //
//          Seoul 151-744, Korea                                           //
//                                                                         //
//          E-mail:  cmp@aces.snu.ac.kr                                    //
//                                                                         //
//-------------------------------------------------------------------------//

//-------------------------------------------------------------------------//
// Authors: Sangmin Seo, Jungwon Kim, Jun Lee, Jeongho Nah, Gangwon Jo,    //
//          and Jaejin Lee                                                 //
//-------------------------------------------------------------------------//

#include "header.h"
#include "mpinpb.h"

//---------------------------------------------------------------------
// block-diagonal matrix-vector multiplication                       
//---------------------------------------------------------------------
void pinvr()
{
  int i, j, k;
  double r1, r2, r3, r4, r5, t1, t2;
  int scount1[total_nodes], rcount1[total_nodes], displ1[total_nodes];
  //double trhs[KMAX][JMAXP+1][IMAXP+1][5]={0.0};
  int kstart, kend,kgap,krem;

  kgap= grid_points[2]/ total_nodes;
  krem = grid_points[2]% total_nodes;
  if(krem == 0){
	kstart=node*kgap;
	kend= kstart+kgap;}
  else {
	  if(node+1 <= krem) 
	  {
		kstart = node*(kgap+1); //1 , 4
		kend= kstart+kgap+1; // 4, 7,
	  }
	  else
	  {
		kstart= krem*(kgap+1)+(node - krem)*kgap; //	7 , 9
		kend = kstart +kgap; // 9, 11
	  }
	}

  for(i=0; i< total_nodes; i++)
  {
    if(i+1 <= krem)
    {
	rcount1[i] = (kgap+1)*(JMAXP+1)*(IMAXP+1)*5;
	scount1[i]= rcount1[i];
    }    
    else
    {
	rcount1[i]= kgap*(JMAXP+1)*(IMAXP+1)*5;
	scount1[i]= rcount1[i];
    }
    if(i == 0)
	displ1[i] = 0;
    else
	displ1[i] = displ1[i-1]+rcount1[i-1]; 	
  }
  if(node == root)
	kstart = kstart+1;
  if(node == total_nodes -1)
	kend = kend -1;	

  if (timeron) timer_start(t_pinvr);
  //#pragma omp parallel for default(shared) private(i,j,k,r1,r2,r3,r4,r5,t1,t2)
  for (k = kstart; k < kend; k++) {
    for (j = 1; j <= ny2; j++) {
      for (i = 1; i <= nx2; i++) {
        r1 = rhs[k][j][i][0];
        r2 = rhs[k][j][i][1];
        r3 = rhs[k][j][i][2];
        r4 = rhs[k][j][i][3];
        r5 = rhs[k][j][i][4];

        t1 = bt * r1;
        t2 = 0.5 * ( r4 + r5 );

        rhs[k][j][i][0] =  bt * ( r4 - r5 );
        rhs[k][j][i][1] = -r3;
        rhs[k][j][i][2] =  r2;
        rhs[k][j][i][3] = -t1 + t2;
        rhs[k][j][i][4] =  t1 + t2;
      }
    }
  }
  if (timeron) timer_stop(t_pinvr);
  
  MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&rhs[0][0][0][0],rcount1,displ1, MPI_DOUBLE,MPI_COMM_WORLD); 
  

}

